Effects of disorder on lattice Ginzburg-Landau model of 
d-wave superconductors and superfluids 
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We study the effects of quenched disorder on the two-dimensional d-wave superconductors (SC's) 
at zero temperature by Monte-Carlo simulations. The model is defined on the three-dimesional 
(3D) lattice and the SC pair field is put on each spatial link as motivated in the resonating-valence- 
bond theory of the high-Tc SC's. For the nonrandom case, the model exhibits a second-order phase 
transition to a SC state as density of charge carriers is increased. It belongs to the universality 
class different from that of the 3D XY model. Quenched disorders (impurities) are introduced both 
in the hopping amplitude and the plaquette term of pair fields. Then the second-order transition 
disappears at a critical concentration of quenched disorder, pc — 15%. Implication of the results to 
cold atomic systems in optical lattices is also discussed. 

PACS numbers: 74.81.-g, 11. 15. Ha, 74.25.Dw 



Effects of disorder on phase structures and phase tran- 
sitions have been studied for various systems. In par- 
ticular, the high-Tc materials are microscopically highly 
nonuniform and it is suggested that there exists a spin- 
glass like phase near the phase transition point of super- 
conductivity (SC) at low temperatures (T')[l[. Further- 
more, existence of a Bose glass was recently suggested 
in the Mott-insulating phase of cold atomic systems in 
random potentials [3]. In the present paper, being mo- 
tivated in part by these observations, we shall study ef- 
fects of quenched disorders on SC's by using a lattice 
Ginzburg-Landau (GL) model of unconventional d-wave 
SC's. The model in its pure case was introduced as the 
GL theory of the resonating-valence-bond (RVB) field 
for the t-J model fsll, and the case of including interac- 
tion with the electromagnetic (EM) field has been inves- 
tigated recently [1]. We expect that this model also de- 
scribes superfluid phase of fermionic atoms in cold atomic 
systems in optical lattices. There, the effects of disorders 
can be investigated well under control. 

We are interested in quantum SC phase transition of 
the two-dimensional (2D) model at T = 0. The model in 
path-integral representation is described in terms of the 
following RVB- type Cooper pair field Uxj- 

Uxj ~ {■4;x+j,'\i^xi - ■ilJx+j,iipx^), (1) 

where x{xo, xi, X2) is the site of the 3D cubic lattice 
of size V = with periodic boundary condition and 
j — 1,2 denotes the spatial direction and also the unit 
vector in j-th direction. Tp^a is the electron annihilation 
operator at site x and spin a I- In the d^2_y2-wave 
SC, the Cooper pair amplitudes in a; = xi and y — X2 
have the opposite signatures as {UX1UI2) < 0. 

Below we shall neglect the effects of the EM field and 
focus on the effects of quenched disorders. The action of 
the clean system is then given as follows: 

C2 Ux2Ux+2,lUx+l,2Uxl 



+C3 {Ux+l,2Uxl + Ux+2,lUx+l,2 

+ Ux2Ux+2,l + UxlUx2) 
+C4 {Ux+2aUx1 + Ux+l,2Ux2) 

+C5 iUx+0,lUxl + Ux+0,2Ux2) + C.C. 



(2) 



where we consider the London limit and set Uxj a U(l) 
variable, Uxj = expiiOxj), Oxj G [— 7r,7r). The overall 
factor g plays a role of 1/fi. and controls quantum fluctu- 
ations. When we fix c^, 5 can be taken as an increasing 
function of the carrier concentration 5. (In the t-J model 
6 is the hole density.) The C2 terms controls fluctuation of 
flux of Uxj^s around each spatial plaquette. The C3 and 
C4 terms represent the spatial hopping of Uxj, whereas 
the C5 term describes the hopping in the imaginary-time 
(0-th) direction. The partition function Z is given by 
Z = J[dU]cxp{A), [dU] =ll^_jdexj/{27r). We consider 

the parameter region C3 < to expect {UxiUl.2) < 0, 
although ^(cs) — Z{—cs) because of the change of vari- 
ables Uxi —Uxi- The action A is related to the action 
^Higgs of the U(l) Higgs gauge theory that is obtained 
from Eq.(I2]) by the replacement Uxj (j)x+jUxj(j)x, where 
(px = exp(i(^2;) is the U(l) Higgs field. ^Higgs is invari- 
ant under time-independent local gauge transformation 
^ 'Px + Oxj ^ Oxj + Xx+j - K where Xx(xi,x2) is 
an a;o-independent function. Actually, A is viewed as the 
gauge-fixed version of ^Higgs in the unitary gauge 4>x = ^■ 
Queched disorder is introduced in the system ^ by 
replacing the coefhcinets 02,3,4 as spatial-plaquette de- 
pendent ones. First, we consider the 2D spatial plane 
with fixed xq, say xq — 0. Among plaquettes in 
the plane we choose p x plaquettes randomly as ones 
at which impurities reside. We call it a sample. We 
consider that the configuration of "wrong" plaquettes 
is a;o-independent because the location of impurities are 
fixed along the imaginary time. Then we reverse the 
values of C2 3,4 for interaction terms contained in these 
plaquettes p. Thus the new plaquette-dependent coefh- 
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FIG. 1: Specific lieat C for p = as a function of g with 
C2 = — C3 = C4 = C5 = 1. Error bars represent MC errors. 
System-size dependence of its peak supports existence of a 
second-order pliase transition. 
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FIG. 2: Correlation function G(r) for various values of g with 
C2,4,5 = — C3 = 1. Critical coupling is estimated gc = 0.154 
for L = 20. Data show that for g > gc there exists SC LRO, 
whereas g < gc no LRO. 



cients 3 4 are given by 



2,3,4 



C2,3,4 with probability 1 - 
— C2,3,4 with probability p, 



(3) 



for each spatial plaquette. Please note — c^. Then the 
partition function Zp and the free energy Fp per site of 
one sample are given by 



[dU] exp(^p) = exp(-yFp 



(4) 



where Ap is obtained from Eq.(l2|) by replacing c; by 
. To obtain an ensemble average (O) of an observ- 
able 0{{Uxj}) in the disordered system, we first prepare 
Ns samples and calculate the quantum-mechanical aver- 
age {0)s for the s-th sample (s — 1, • • • , Ng). Then we 
average it over samples, 



(O) 



1 ^= r 

— ^(O)., (O), - J [dU]OeMAp. 



(5) 



For the MC simulations, we used the standard 
Metropolis algorithm 7] . The typical statistics used was 
10^ sweeps per block and the average and MC errors 
were estimated over 10 blocks for each sample. Then we 
take quenched averages over Ng = 30 ^ 50 samples. We 
estimated standard deviation of physical quantities like 
"specific heat" over samples (we call it sample error) as 

a function of Ns, which becomes stable for Ng ^ 30. 

Let us first consider the nonrandom case, p = 0. We 
studied the phase structure by calculating the "internal 
energy" per site E = —{A)/V and the "specific heat" per 
site C = {{A-{A))'^)/V. Typical behavior of C is shown 
in Fig[Tl which indicates a second-order phase transition 
ai g = Qc — 0.154. This transition has been predicted in 
Ref. Pi and also observed in the presence of the EM gauge 
interactions 4] . In order to verify that it is a transition 
to a SC phase, we measured the correlation function of 



Uxj, the order parameter field of SC, 



(6) 



We present G{r) in FigO It is obvious that there ex- 
ists SC long-range order (LRO) for g > g^ whereas no 
LRO for g < 5c [8j. Therefore the observed transition is 
nothing but the SC phase transition. 

The critical exponents for 02,4,5 = —C3 = 1 were es- 
timated by the finite-size scaling analysis for C. We 
obtained i' — 1.5, a — 0.285 and the critical coupling 
goo = 0.153. When 0-2 — 03 — 0, the system © reduces 
to a set of decoupled 2D XY spin models [Uxj plays the 
role of an XY spin], each of which has the Kosterlitz- 
Thouless transition. The C2 and C3 terms couple these 
2D XY spins. From the above values of critical expo- 
nents, we judge that the present phase transition does 
not belong to the universality class of the 3D XY model 
having v = 0.6721 (13) [9]. 

Let us next turn to the random case, p > 0. At first, 
we consider the gc2 — gc' plane where c' = — C3 = 04 = 05, 
and search for the location of the peak of C. In FiglJl 
we present the peak location of C for given disorder con- 
centration p = 0.0,0.10,0.20. For p — 0.0, the peak line 
expresses the second-order transition as we saw in FigsU] 
and [21 We see that the region of the normal (non SC) 
state increases as p increases, although we need to check 
whether the location of this peak expresses genuine SC 
transition for p > 0. Below we examine it by focusing on 
the specific case C2 = — C3 = C4 = C5 = 1 with the varying 
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FIG. 3: Location of the peak of C with L = 12 and iV^ = 30 
in the gc2 — g x (— C3 — C4 = C5) plane for p = 0,0.10,0.20. 
The region of disorder "phase" increases as p increases. 
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FIG. 4: Specific heat C for L = 16 with A^, = 30 vs g for 
p = 0.0, ■ ■ ■ , 0.20. Error bars here and below represent sample 
errors. The peak around the SC phase transition is getting 
weaker as p is increased. Location of the peak moves to larger 

g- 



parameter g and calculate E, C, G{r), etc. Density of 
quenched disorder that we studied is p = 0.05,0.10,0.15 
and 0.20. 

We first present the result of C in Figd) The signal 
of the second-oder phase transition at p = is getting 
weaker as p is increased, and also the location of the 
peak moves to larger g. 

Next, in Fig[51 we present the system-size dependence 
(SSD) of C for p = 0.10, 0.15 and 0.20. We also calcu- 
lated the derivative of C, D{g) = gdC{g)/dg, to identify 
the order of the phase transition. Results are shown in 
Figini From these results, we judge that there exists a 
second-order transition for p < 0.10 whereas it changes 
to a crossover for p > 0.20. The case oip — 0.15 seems to 
have no SSD in C but have certain SSD for D{g). There- 
fore, the transition in the p = 0.15 case might be of third 
order. 

To study if the SC state persists even for p > 0.20, we 
measured the SC correlation G{r) averaged over samples 
with randomly generated c^. In Fig[7] we present the 
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FIG. 6: D{g) = gdC(g)/dg with A^^ = 30 vs 5 for various L. 
(a)p = 0.10; (b)p = 0.15; (c)p = 0.20. The cases (a) and (b) 
show SSD, while the case (c) shows almost no SSD. 



results for p = 0.10, 0.15 and 0.20. It supports that for 
p = 0.10 there exists the SC LRO for g > gc whereas for 
p = 0.20 no LRO for any values of g. This observation 
leads to the conclusion that the SC phase disappears for 
p > 0.20, i.e., there is a multicritical point p^^ near p = 
0.15 in the g — p phase diagram. 

It is interesing to ask what kind of phase is realized 
in the region p > pmc and g > gcr, where ^cr is the 
crossover coupling determined by the peak location of the 
specific heat. Probably, in the case of the doped AF mag- 
nets with strong inhomogeneity, this phase may simply 
correspond to dirty "normal" metallic state. We think 
that the present model given by Eq. ^ also describes the 
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FIG. 5: The system size dependence of C with Na = 30 for 
p = 0.10, 0.15, 0.20. C for p = 0.10 shows SSD, which 
supports a second-order phase transition to SC state. C for 
p = 0.15 has less SSD, and C for p = 0.20 shows no SSD for 
L = 16, 24. 



FIG. 7: SC correlation function G{r) for L = 20 with A^^ = 50. 
(a)p = 0.10; (b)p = 0.15; (c)p = 0.20. (d) G(r = 10) for 
(a,b,c) vs g. SC "transition point" gc is estimated from the 
specific heat C as gc ~ (a)0.168, (6)0.17, (c)0.18. Fig.(d) 
supports that LRO developes for g > gc in (a). 
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FIG. 8: Phase diagram in the p — phase diagram for the 
lattice model of dirty d-wave SC with C2 = — cs = C4 = cs = 1. 
Possible Bose glass phase is suggested. There should be a 
multicritical point near the line p — 0.15. 



superfluid d-wave RVB state of fermionic atoms in 2D 
optical lattice, which was recently proposed by sereval 
authors p^. Random disorders can be introduced into 
the systems by a laser speckle or by an incommensurate 
bichromatic potentail. In that experimental setup, there 
is an interesting possibility that a Bose glass phase, which 
is an analogue of the spin glass phase, is realized there. 
For ultracold strongly interacting ^^Rb bosons, a Bose 
glass phase is suggested by experiments [1]. In the ex- 
periments of cold atom systems in optical lattices, the 
Bose glass phase has no long-range coherence but exci- 
tations are gapless. In FiglSl an expected phase diagram 
of the present lattice model of dirty d-wave SC and su- 
perfluidity (SF) is shown. We think that existence of the 
Bose glass phase is examined theoretically by the stan- 
dard analytical and numerical methods applied for spin 
glass [nl [13. 
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